Uncertainties in Pion and Kaon Fragmentation Functions 
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We present a detailed assessment of uncertainties in parton-to-pion and parton-to-kaon fragmen- 
tation functions obtained in recent global QCD analyses of single-inclusive hadron production data 
at next-to-leading order accuracy. We use the robust Lagrange multiplier approach for determining 
uncertainties to validate the applicability of the simpler but approximate Hessian method. Extensive 
comparisons of the results obtained within both methods are presented for the individual parton- 
to-pion and kaon fragmentation functions. We provide Hessian eigenvector sets of pion and kaon 
fragmentation functions that allow one to easily propagate their uncertainties to any observable. 
Various applications of these sets are presented for pion and kaon production in electron-positron 
annihilation, lepton-nucleon scattering, and proton-proton collisions. 
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I. INTRODUCTION AND MOTIVATION 

Fragmentation functions (FFs) contain vital non- 
perturbative information required to describe a great va- 
riety of hard scattering processes with a given identified 
hadron in the final-state within the framework of per- 
turbative Quantum Chromodynamics (pQCD). The en- 
hanced sensitivity of such less inclusive probes to crucial 
features of the partonic structure of hadrons, including 
its flavor content, dependence on spin [J Hj], or possible 
nuclear modifications 0, H|, combined with ever grow- 
ing experimental precision, opens up novel opportuni- 
ties for phenomenological QCD studies with an accu- 
racy hitherto reserved to fully inclusive measurements. 
In this context, the availability of reliable sets of parton- 
to-hadron FFs, as well as accurate estimates of their un- 
certainties, is of the utmost relevance. 

The extraction of FFs from data, originally restricted 
to single-inclusive electron-positron annihilation (SI A) 
has evolved in recent years into truly global QCD 
analyses encompassing data obtained in hadron-hadron 
collisions as well as hadron multiplicities in semi- 

inclusive deep-inelastic scattering (SIDIS) @. Only the 
combined analysis of the large body of existing data taken 
in different processes and for different hadron species 
allows for a consistent determination of all aspects of 
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hadronization such as fully charge and flavor separated 
parton-to-hadron FFs. A global fit also constitutes an 
explicit check of the assumed underlying factorizability 
and universality of FFs in hard processes, which is the 
foundation of the predictive power of pQCD. 

As is the case with parton distribution functions 
(PDFs) [lH, the assessment of uncertainties for FFs ob- 
tained in global QCD analyses is far from being straight- 
forward. Such studies incorporate data from many differ- 
ent experiments with diverse characteristics and errors. 
The data are in turn confronted with theoretical esti- 
mates whose inherent uncertainties are notoriously diffi- 
cult to quantify, because, in addition to the truncation 
of the perturbative expansion at a given order of pQCD, 
there are also unavoidable approximations and assump- 
tions involved. The latter comprise the choice of the 
functional form used to parametrize FFs 13] (or PDFs 
14] ) at some low input scale Q — 1 GeV, the value of 
the strong coupling at some reference scale, the selection 
of and cuts applied to the data sets used in the fit, or the 
treatment of heavy flavors in the scale evolution and the 
calculation of observables. 

In recent years, significant progress has been made 
in putting forward, exploring, validating, and compar- 
ing different strategies to estimate uncertainties in global 
QCD analyses of PDFs. Among the various approaches, 
the robust Lagrange Multiplier (LM) technique [ll| di- 
rectly relates the variation of the parameters determined 
in the fit or, more generally, of any observable computed 
with them, to the variation of the x 2 function that quan- 
tifies the goodness of the fit to data. While this avoids 
any approximations or assumptions about the depen- 
dence of the \ 2 hypersurface on the parameters used to 
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describe the PDFs or how to propagate uncertainties to 
a given observable, it requires an extensive amount of x 2 
minimizations involving all data sets included in the fit. 
The computationally less demanding "Improved (itera- 
tive) Hessian" (IH) approach [16| assumes a quadratic 
behavior of the x 2 hypersurface on the parameter dis- 
placements and represents the x 2 increment from its min- 
imum value in terms of combinations of the fit parame- 
ters that maximize the variation. Within this eigenvec- 
tor representation of the Hessian matrix, and assuming 
a linear propagation of errors, the uncertainty of any ob- 
servable can be straightforwardly estimated from a set of 
pre-calculated fits corresponding to fixed displacements 
along the eigenvector directions. Another approach is 
based on analyzing a large amount o f rep licas of the orig- 
inal data sets with neural networks 17]. While this has 
the drawback of defining the central fit only as the sta- 
tistical average of 0(100) PDF fits, it is largely free of 
the bias from assuming a certain functional form for the 
PDFs at scale Q a . 

Corresponding studies of uncertainties in parton-to- 
hadron FFs are scarce despite the large body of available 
data. An analysis based solely on SIA data, which do 
not allow for a full charge and flavor separation of FFs, 
was presented in Ref. |8[ using the Hessian method but 
without validating its applicability. The so far most com- 
prehensive analysis based on data from SIA, SIDIS, and 
hadronic collisions was presented in Refs. || hence- 
forth referred to as DSS FFs. However, uncertainties 
of the extracted FFs were only assessed qualitatively for 
certain truncated moments within the LM approach. 

In the following, we will extend the DSS analysis by 
performing a detailed assessment of the uncertainties in 
pion and kaon FFs based on the IH method. We closely 
follow the DSS framework Q and adopt the same func- 
tional form to parametrize the FFs at the initial scale 
Qo = 1 GeV and the same selection of data sets. In gen- 
eral, we find good agreement with the results obtained 
with the robust LM technique but notice that for the 
same nominal tolerance criterion, i.e., increase in x 2 , the 
IH method typically leads to somewhat smaller uncer- 
tainty estimates. We devise a recipe to account for these 
small differences and provide sets of Hessian eigenvectors 
FFs to facilitate the propagation of uncertainties to ar- 
bitrary observables. We believe that such an analysis is 
particularly timely and useful in view of the wealth of up- 
coming or new precise data on identified hadron yields in 
SIA from B factories QJ], SIDIS and hadron-hadron 
collisions both at RHIC [H and the LHC [H|. The re- 
sults of our analysis will help to quantify the impact of 
these data sets, identify possible tensions with the DSS 
analysis, and will serve as the baseline result for an an- 
ticipated update of the DSS sets of pion and kaon FFs. 

The remainder of the paper is organized as follows: 
in the next Section we briefly recall the main aspects of 
the DSS analysis and the IH method, study in detail the 
uncertainties of parton-to-pion FFs in IH approach, and 
compare to the results obtained with the LM technique. 



Hessian eigenvector sets of pion FF are provided and ap- 
plied to calculations of pion yields in SIA, SIDIS, and 
hadron-hadron collisions. Section III is devoted to a sim- 
ilar study for kaon FFs. We briefly summarize the main 
results in Sec. IV. 



II. PION FRAGMENTATION FUNCTIONS 
A. Preliminaries 

Since both the framework and methodology for the ex- 
traction of FFs in global QCD analyses at next-to-leading 
order (NLO) accuracy as well as the implementation of 
the IH and LM techniques have already been explained 
in quite some depth in the literature d, [IH, HH , we will 
only briefly recall the main concepts and results relevant 
for our studies. 

As was mentioned in the Introduction, we choose the 
DSS NLO analysis of pion and kaon FF [9( and their 
error estimates based on the LM method as our baseline 
fit. We adopt the same selection of data sets used in the 
DSS fit and the same flexible functional form 

T~)H f n v = A^(f-Z)ft[f+ 7Z (I-Z)^ 

1 K ' W0) B[2 + a i ,pi + l]+>Y i B[2 + ai,pi + 6 i + l]' 

to parametrize the hadronization of a parton i into a 
hadron H at the initial scale Qo = 1 GeV. z denotes the 
fraction of the parton's momentum taken by the observed 
hadron, and B[a, b] represents the Euler Beta- function. 
The Ni in ([T]) are normalized such to represent the con- 
tribution of Df 1 to the momentum sum rule 

V / dzzD?(z,Q 2 ) = l. (2) 

H J ° 

Since presently available data do not constrain all free 
parameters in (JT]) for each parton i equally well, certain 
relations upon the individual FFs had to be imposed in 
the DSS analysis @ without jeopardizing the quality of 
the fit. Apart from assuming isospin symmetry for the 
unfavored sea quark FFs, i.e., Df = DJ , the total 
w-quark £>jtj and d-quark Df^_ a FFs are only allowed 
to differ in normalization TV, and the strange quark FFs 
Df = Df are related to the sea quark FFs through 
another normalization factor TV'. The corresponding FFs 
for negatively charged pions are obtained by charge con- 
jugation invariance and those for neutral pions by assum- 
ing Df = [Df + + Df }/2. 

The remaining 23 free parameters {a{\ describing the 
DSS FFs for quarks and gluons into positively charged 
pions, Df , are then determined by a standard x 2 mini- 
mization. Since the full error correlation matrices are not 
available for most of the data sets used in the fit, statisti- 
cal and systematical errors are usually added in quadra- 
ture ■ As in the DSS analysis, we allow each data 
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set to float within the quoted experimental normalization 
uncertainty by introducing a set of 7 extra parameters to 
the fit. 

It suffices to say that we fully reproduce the set of 
optimum parameters {a?} of the DSS analysis [9|, corre- 
sponding to the minimum in the x 2 profile. To explore 
their uncertainties with the IH approach [16[ we express 
the Hessian matrix 



y 2 dyidy-j 



(3) 



where the derivates are taken at the minimum, in terms of 
its A par eigenvectors . The displacements j/j = a; — a" 
in Eq. (|3]) and the increase in x 2 

A X 2 = X 2 ({a t }) - X 2 o({a°}) = £ H ijyiVj (4) 



are then replaced by a new set of parameters {zi} defined 
by M 
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(5) 



The factors sj oc y/l/sj are used to rescale the {z{\ such 
that the distance from the x 2 minimum is simply given 
by 



Ax 2 = £. 



(6) 



The x 2 function changes rapidly for directions corre- 
sponding to large eigenvalues £k of the Hessian matrix 
while small eigenvalues belong to directions where the fit 
parameters are only weakly constrained. 

What makes the eigenvector representation {zt} par- 
ticularly useful and convenient is the possibility to con- 
struct 2 A par basis sets Sj^ of FFs which greatly facilitate 
the propagation of their uncertainties to arbitrary ob- 
servables O. An estimate of the error AO away from its 
best fit estimate O(S ) is obtained by computing 16 | 



AO = - 

2 



k=l 



1/2 



(7) 



The eigenvector sets 



are defined by choosing the 



amount T = y A\ 2 still tolerated for an acceptable 
global fit and correspond to positive and negative dis- 
placements by T along each of the eigenvector directions 



±TSi, 



(8) 



We will extensively use the obtained basis sets and 
Eq. (J?} in the remainder of the paper to propagate un- 
certainties of FFs to several observables. 
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eigenvector directions {z.} 

FIG. 1: Correlations between the fit parameters {a,i} and the 
eigenvector directions {zi}. The larger the box size the larger 
the overlap, see text. 



B. Results 

First, before making use of the eigenvector sets, we 
need to validate the Hessian method as a reliable tool 
to determine uncertainties of pion FFs. Figure [1] shows 
the overlap of each of the original fit parameters {a^} in 
Eq. ((T|) with the eigenvector directions {z{\ introduced in 
Sec. Ill Al The larger the box size the larger the contribu- 
tion of a certain eigenvector direction to a fit parameter 
a*. The {z;} in Fig. [1] are ordered in terms of the size 
of the corresponding eigenvalues of the Hessian matrix: 
21 corresponds to the largest eigenvalue, i.e., the direc- 
tion in parameter space where % 2 changes most rapidly, 
whereas 230 is only very weakly constrained by data. One 
can see that in most cases there is a fairly strong correla- 
tion between a given original fit parameter a.; and a single 
eigenvector direction Zj. The opposite case, when several 
fit parameters are strongly correlated with an eigenvector 
direction, implies that those fit parameters are mutually 
correlated cannot be constrained independently. 

As can be inferred from Fig. [I] the relative normal- 
izations applied to data sets in the fit, {A^aleph, Asld, 
Arermes, A STA r, A T pc, A PH enix, A pal}, are among 
the best constrained parameters as they are typically 
linked to only a very small number of eigenvector direc- 
tions, all corresponding to large eigenvalues. The range 
of variation of these parameters is governed by the nor- 
malization uncertainty quoted by each of the experiments 
and their main role is to ease possible tensions among the 
data sets in the fit. 

In a second group of comparatively well constrained 
parameters one can find the normalization factors Aj, 
see Eq. ([1}, of the different FFs for positively charged 
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FIG. 2: [color online] Examples of deviations from the ex- 
pected parabolic behavior A% 2 = T 2 for selected eigenvector 
directions {zt}, see text. 

pions Df + , starting from the dominant (or favored) one, 
N u +u, followed by N g for the gluon FF, N b+ i, N a , and 
N c+S . The Ni represent the second moments of the FFs 
for flavor i into hadron H, 

D?{Q 2 ) = f zD?{z,Q 2 )dz, (9) 
Jo 

entering the sum rule ([2]). One notices that the different 
Ni are strongly correlated with each other. 

The parameters ai and Pi, controlling the main fea- 
tures of the z dependence in Eq. ([1} , along with N and 
N' , related to certain flavor symmetry relations among 
the FFs, fall into the next category, all showing fairly 
large correlations. Finally, the parameters associated 
with more subtle details of the z dependence of the FFs, 
7i and Si, are the least well constrained ones in the fit 
and are mainly correlated with eigenvector directions be- 
longing to smaller eigenvalues. 

Even though the correlations shown in Fig. Q] make ex- 
plicit the hierarchy of the fit parameters with respect to 
the level of how well they are constrained and to what 
extent correlations among them can be found, it does not 
necessarily indicate if x 2 exhibits the assumed quadratic 
behavior on the parameters away from the best fit. To 
explore this further, Fig. [5] illustrates the deviations of 
the x 2 function from the expected quadratic dependence 
for selected, representative eigenvector directions {zi}. 
Here, we vary only one of the parameters z% at a time 
such that a given change of A% 2 = T 2 is produced. Of 
course, since each Zi has in principle overlap with all fit 
parameters {ai}, the latter all vary in this procedure. 
For a truly quadratic behavior near the minimum, as is 
the underlying assumption in the Hessian approach, the 
quantity T 2 — A\ 2 vanishes, where Ax 2 is the change 



in x 2 induced by the variation of the parameter Zj. Any 
deviation from zero will signal a departure from the ex- 
pected parabolic dependence. 

As can be seen, for all the parameters Zi shown in 
Fig. [2] the choice T = 1, i.e., A% 2 = 1, works very well, 
leading only to fairly small deviations from zero. This 
implies that Hessian method is reliable for Ax 2 = 1, and 
our eigenvector sets are expected to produce faithful 
uncertainty estimates close to those obtained with the 
robust LM approach. For some of the eigenvector direc- 
tions Zi, those strongly constrained by data and corre- 
sponding to large eigenvalues of the Hessian matrix, the 
quadratic behavior persists even far away from the x 2 
minimum, i.e., up to large values of T. Examples are Z3 
and zq, which are mainly correlated with the normaliza- 
tions N u+a and N c+S of the total up and charm quark 
FFs, respectively. Some eigenvector directions, like, for 
instance, Zxg, according to Fig. [T] mainly related to the 
small z exponent a a of the anti-up-quark FF, and Zig, 
controlling the large z behavior of the total charm quark 
FF through /3 c +g, do show a more pronounced depar- 
ture from the ideal parabolic behavior starting already 
at about T 2 ~ 5. Not surprisingly, deviations are most 
pronounced for poorly constrained parameters such as 
Z25 and Z30 mainly correlated with the dependence of to- 
tal up quark and gluon FFs at intermediate values of z 
as described by 6 u + a and j g , respectively. 

The deviations from the assumed quadratic behavior 
observed above are also reflected in the actual x 2 P ro ~ 
files for the fit parameters mainly associated with the 
eigenvector directions z, shown in Fig. [5J They are read- 
ily computed within the LM technique and displayed in 
Fig. [3] Clearly, while for some of the parameters the 
profiles are reasonably smooth and parabolic, as is as- 
sumed in the Hessian method, in general they are not. 
The profiles can exhibit asymmetric shapes, almost flat 
regions, and perhaps even multiple minima. It is worth 
pointing out that in case of the DSS fit [9j none of these 
features is related to a lack of flexibility in the chosen 
functional form, Eq. ([T]). For instance, tensions among 
different data sets cause the asymmetry in the x 2 profile 
for a-u, and the flat part in the profile for 7 9 is caused 
by insufficient constraints from data in some kinematic 
regions. 

At this point it is important to check whether the 
IH approach yields meaningful uncertainty estimates for 
observables, despite the deviations from the assumed 
quadratic dependence of x 2 on the parameters, illus- 
trated in Figs. [5] and [3J In the DSS analysis [9j, truncated 
second moments of the FFs were identified as representa- 
tive indicators of the typical uncertainties in the fit, and 
their variations were studied within the LM framework. 
They are defined in complete analogy to Eq. $2$ by in- 
troducing a lower cut-off z m i n to the range of integration: 

77fOmin,Q 2 )= / zD?( Z ,Q 2 )dz. (10) 

To avoid the kinematic region of small z where mass 
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FIG. 3: The \ 2 profiles obtained with the LM approach for 
the six fit parameters discussed in Fig. [2] The horizontal lines 
indicate an increase A% 2 by one unit or by 2% of the total \ 2 
as adopted in the DSS analysis [ij; see text. 

effects, neglected in FFs, become relevant, z m i n = 0.2 
was chosen in [9]. We note that the picture arising 
from the DSS analysis of the truncated second moments 
f]f (^mim Q 2 ) [2( was recently shown to agree well with ex- 
plicit calculations of uncertainties in hadron production 
cross sections at LHC kinematics based on the robust LM 
method [H- 

In Figure @] we show the x 2 -P ron l es (solid lines) for 
the truncated second moments rj^ + of the DSS pion FFs 
at NLO accuracy at Q = 5 GeV and for selected parton 
flavors i = {u + u, g,c + c} as obtained with the LM 
technique. The moments are normalized to the value 770 
they take for the optimum fit to data, characterized by 
the set of parameters {a°}. From the width of the curves 
at a given value of A\ 2 , one can read off the uncertainties 
on the various rj/rjo, which can be in general asymmetric 
with respect to the best fit, i.e., ?y/?/o = 1- 

In the DSS analysis of pion and kaon FFs |9| , a typi- 
cal tolerance of A\ 2 /x 2 — 2% was regarded as a faithful 
measure of uncertainties, which amounts to an increase 
in the total x 2 of the fit to pion data by about 17 units 



as is indicated by the horizontal lines in Fig. [4] labeled as 
"LM 2%" . It is important to recall that this tolerance cri- 
terion is not derived from any theoretical argument but 
is an empirical estimate by requiring that all data sets 
included in the global fit are still adequately described. 
Such kind of choices at various levels of sophistication are 
typically made also in global QCD analyses of polarized 
and unpolarized parton densities, see, e.g., [l], HJ. The 
naive criterion A\ 2 — 1 is usually regarded as too small 
to reliably account for PDF or FF uncertainties due to 
the complex nature of global fits, the different charac- 
teristics of the data sets, and various, often unaccounted 
sources of non-Gaussian theoretical errors. As can be 
inferred from Fig. HJ variations of the truncated second 
moments are typically found to be of the order of 2 -f- 3%, 
5%, and around 10% for the u + u, gluon, and c + c to 
pion FFs. 

Also in Fig. [H the horizontal lines labeled as "IH 2%" 
and "IH 5%" represent the same uncertainties as ob- 
tained with the LM method above but now are estimated 
within the IH framework for two different tolerances (2 
and 5 percent increase in the total \ 2 of the fit, respec- 
tively). Clearly, the 2% criterion systematically under- 
estimates the uncertainties found in the robust LM ap- 
proach. Not surprisingly, the differences have their origin 
in the assumed quadratic behavior of the % 2 profile away 
from its minimum, which, in general, is not fully ade- 
quate as we have already demonstrated above in Fig. [2l 

Since the discrepancies between the LM and IH es- 
timates based on Ax 2 /x 2 = 2% are not large enough 
to completely invalidate the Hessian method, we devise 
a simple recipe to remedy its shortcomings in practi- 
cal applications based on pion FFs. We observe, that 
the IH method reproduces the LM results obtained for 
^X 2 /x 2 — 2% much better if one allows for a larger tol- 
erance of 5%, as is indicated by the horizontal lines for 
"IH 5%" in Fig. Q] Hence, we propose to use preferably 
the IH method with Ax 2 /x 2 = 5% to closely match the 
uncertainty estimates for pion FFs advocated in the DSS 
analysis [9|] and provide relevant Hessian eigenvector ba- 
sis sets S± for both a A X 2 /x 2 of 2% and 5% 

As an application and also to further verify the use- 
fulness of the IH method, we calculate uncertainty es- 
timates for various representative data sets used in the 
DSS analysis of pion FFs. To propagate the uncertain- 
ties of the FFs to pion production cross sections we use 
Eq. (O throughout. Figure [5] illustrates the agreement 
between preliminary data for charged pion multiplicities 
[24j in SIDIS and the corresponding results of the DSS fit 
at NLO accuracy [9( in various bins of z. The results are 
largely independent of the choice of PDFs. The multiplic- 
ity data are instrumental in providing flavor and charge 
separated quark-to-pion fragmentation functions. We no- 
tice that the differences between uncertainties obtained 
with the 2% and 5% Hessian sets are less pronounced 
than in Fig. [U suggesting that the 2% variations already 
account for most of the error in the kinematic region 
probed by the multiplicity data. 
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FIG. 5: left: comparison of the NLO DSS results for charged pion multiplicities in SIDIS with preliminary data [24[ as a 
function of Q 2 in various bins of z. right: "(data-theory)/theory" for the DSS results; open and full circles denote 7r + and tt~ 
multiplicities, respectively. For all panels, the light and dark shaded bands indicate uncertainty estimates for A^ 2 /x 2 = 2% 
and 5%, respectively, based on our Hessian eigenvector sets and using Eq. (O. 



Similar comparisons for single-inclusive pion produc- 
tion in SIA and pp collisions at BNL-RHIC are shown in 
Fig. [6] Data are taken from [25| and [26|, respectively. 
Again, the differences between the uncertainty estimates 
obtained with the 2% and 5% Hessian eigenvector sets 
are minimal, except for regions sensitive to large momen- 
tum fractions, z > 0.6, where experimental constraints 
become very scarce. In these regions, uncertainties are 
severely underestimated and not trustworthy. In case of 



pp collisions (right-hand-side of Fig. we include for 
comparison an estimate of the theoretical ambiguity due 
to the choice of the factorization scale in the NLO calcu- 
lation [27) which is much more significant than errors on 
FFs (or PDFs as was shown in Ref. [IH). 

A common feature of the results shown in Figs. [5] and 
|5] is that the estimates of the relative uncertainties due 
to FFs remain almost constant in a wide range of hadron 
momentum fractions z and energy scale fx set by the pho- 
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FIG. 6: Uncertainty estimates for single-inclusive pion production in SIA [2^| (left) and in pp collisions at 200 GeV and mid- 
rapidity [2(| (right). The light and dark shaded bands indicate the results for A\ 2 /x 2 = 2% and 5%, respectively, based on 
our Hessian eigenvector sets and using Eq. For pp — > ir°X, we also show the theoretical scale ambiguity (outermost bands). 





FIG. 7: Uncertainty bands for the DSS pion FFs estimated with the IH and LM methods at Q 2 = 10 GeV 2 and Q 2 = M 2 Z . 



ton virtuality Q in case of SIA or SIDIS or the transverse 
momentum pr inpp collisions. As expected, the obtained 
uncertainty bands cover and reflect the typical range of 
the statistical errors of the fitted data relative to the best 
fit prediction. We note that results very similar to those 



shown in Figs. [5] and [6] have been obtained for other ob- 
servables depending on pion FFs. 

Finally, we address the uncertainty estimates on the 
individual parton-to-pion FFs. Figure [7] shows the NLO 
DSS zDf (z,Q 2 ) for i — u + u, u, s + s, g, c + c, 
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and b + b for two different scaies Q 2 = 10 GeV 2 and 
Q 2 = M| along with our estimates of their uncertain- 
ties using both the IH and LM method. For better visi- 
bility and to facilitate comparisons between the different 
FFs, the lower panels for each flavor show the relative 
uncertainties which are typically of the order of 5% for 
the favored quark combination u + u and around 10% for 
unfavored quark-to-pion and gluon FFs, fairly indepen- 
dent of the scale Q. Uncertainties increase significantly 
for large momentum fractions, z > 0.6, where current 
experimental constraints are insufficient. 

In the DSS analysis [9j, the fragmentation of charm and 
bottom quarks into charged pions is included discontin- 
uously as massless partons in the scale evolution above 
their MS "thresholds", Q — m c ^, with m Cj h denoting the 
mass of the charm and bottom quark, respectively. This 
simplified treatment of heavy flavors is at variance with 
current PDF fits [HI, [HI j where very elaborate schemes 
have been developed to properly include mass effects near 
threshold and to resum potentially large logarithms for 
Q 2 ^> m 2 b . Only flavor-tagged SI A data constrain the 
charm and bottom FFs in the DSS analysis, and their 
uncertainties, albeit reasonably small, should be taken 
with a grain of salt. 

Concerning the comparison between the two methods 
to obtain our uncertainty estimates, we find good agree- 
ment if the IH eigenvector sets with the 5% tolerance 
criterion for A\ 2 /x 2 are used. Only for z < 0.3 and un- 
favored, i.e., less well constrained FFs, the LM method 
yields somewhat larger uncertainties compared to the 
corresponding bands obtained with the IH technique as 
can be seen, for instance in the panels showing the u, 
s + s, and the gluon FFs. 

While the robust LM method is straightforwardly ap- 
plied to cross sections or truncated moments of FFs, it 
is not very handy in determining the z dependent uncer- 
tainties. One can make use, for instance, of the profiles 
shown in Fig. [3] for the variations of each fit parameter 
{oj} to compute the spread in the FFs for each parton 
flavor i. In doing so, one neglects, however, correlations 
among the variations of the fit parameters, which can be 
non-negligible as we have demonstrated in Fig. [T] This 
might lead to an overestimate of uncertainties for the z 
dependent FFs. On the other hand, the Hessian method 
is particularly simple, using Eq. ([7]), provided its appli- 
cability has been carefully established and appropriate 
eigenvector sets have been generated [23] ]. This was the 
purpose of our studies for pion FFs presented in this Sec- 
tion. 



III. KAON FRAGMENTATION FUNCTIONS 

A. Preliminaries 

Next, we proceed with a similar error analysis for the 
NLO DSS kaon FFs 9]. Since even the most precise kaon 
production data in SIA exhibit experimental uncertain- 



ties typically at least twice as large as those found for 
pions, one must expect much less, but still reasonably 
well constrained parton-to-kaon FFs. 

To account for the phenomenological expectation that 
the formation of secondary ss pairs should be suppressed 
in the production of, say, a \K + ) = \us), the two favored 
quark combinations D^ +u and D^ +s are fitted indepen- 
dently in the DSS analysis using the functional form ([1}. 
Indeed, in line with that expectation, the DSS fit prefers 
Df +g > D^ +a . Since presently available data do not 
fully constrain all unfavored kaon FFs, is was assumed 
that they all share the same functional form, i.e., 

D r=Dr=Dr=Df. ( n) 

Corresponding FFs into K~ are obtained, as usual, by 
charge conjugation symmetry, which leaves a total of 
24 free fit parameters {at] describing the DSS FFs for 
quarks and gluons into positively charged kaons. Six ad- 
ditional parameters control the relative normalization of 
the data sets in the DSS analysis §. 

B. Results 

As for the pion FFs, we start our analyses by show- 
ing the correlations between the fit parameters {a^} and 
eigenvector directions {z;} in Fig. [8] As before, the {zi} 
are ordered by the size of the corresponding eigenval- 
ues of the Hessian matrix, and correlations can be found 
across the entire parameter space. In most cases, there 
are, however, fairly strong correlations relating a certain 
eigenvector direction with just a one (or only a few) of 
the original fit parameters. Again, the best constrained 
fit parameters are those related to relative normalizations 
applied to the data sets in the fit. In a second group are 
the normalization factors Ni of the different FFs for pos- 
itively charged kaons Df , see Eq. (Q]). Among the least 
well constrained parameters are mainly those associated 
with subtle details of the z dependence of the FFs, i.e., 
7i and Si in Eq. ([1]). 

In the DSS analysis of kaon FFs [9j , uncertainties were 
estimates within the LM method for the truncated mo- 
ments (fit))) and found to be at least twice as large as 
those for pions. Again, an increase of A% 2 /% 2 = 2%, 
corresponding to about 8 units in x 2 , were regarded as a 
faithful estimate of the typical uncertainties. As a rep- 
resentative example, Fig. |H] shows the corresponding \ 2 
profile (solid line) and error estimate (horizontal bar la- 
beled as "LM 2%") for the dominant total strange quark 
FFs into positively charged kaons. As in Fig. 21 the vari- 
ation of r]^ +s is normalized to its value 770 taken for the 
optimum set of parameters {a°} and is found to be of 
order of about 10%. 

When compared to the uncertainties obtained within 
the IH approach, we find that the results corresponding 
to A\ 2 /x 2 — 2% and 5% slightly under- and overesti- 
mate the range of variation of r]ff + g determined with the 
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and b + b for two different 



Ml along with our esti- 
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eigenvector directions {z.} 

FIG. 8: As in Fig.[T]but now showing the correlations between 
the fit parameters {<n} and the eigenvector directions {zi} for 
positively charge kaons. 
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FIG. 9: As in Fig. [4]but now comparing the range of variation 
(indicated by the horizontal bars) of the truncated moment 
r]f of the dominant (s + s)-to-kaon FFs at Q — 5 GeV; see 
text. 



LM method. Qualitatively very similar results can be 
obtained for all the other parton-to-kaon FFs. Again, we 
provide Hessian eigenvector basis sets S k for both toler- 
ance criteria [23j to facilitate estimates of uncertainties 
for arbitrary observables involving kaon FFs. 

Before discussing an application of our Hessian sets, we 
first study the uncertainty bands for the individual NLO 
DSS parton-to-kaon FFs. Figure HU1 shows zDf + (z,Q 2 ) 



for i = u + u, u, g, s + s, c + c, 
scales Q 2 = 10 GeV 2 and Q 2 - z 
mates of their uncertainties using both the IH (5%) and 
LM (2%) method. As for pion FFs, the lower panels for 
each flavor show the relative uncertainties obtained with 
the IH approach. 

Despite the sizable disparity in size, the two favored 
s + s and u + u FFs exhibit very similar relative uncer- 
tainties: close to 20% at intermediate values of z and a 
significant growth for both large and small z. All unfa- 
vored quark FFs, such as , are not well constrained 
by data and show extremely large uncertainties of 50% 
or more. All these features can be traced back to the 
much larger experimental uncertainties for flavor-tagged 
SIA [25| and SIDIS [24| data for kaons compared to those 
for identified pions. 

Rather unexpectedly, the gluon-to-kaon FF shows 
rather small relative uncertainties over most of the z- 
range, comparable to those for the favored quark FFs. 
This feature is found both within the LM and IH method 
and hence cannot be attributed to the approximations 
made in the latter approach. Since the pp data, which 
mainly constrain at z ~ 0.5, are rather sparse and 

less precise than for identified pions, the result most likely 
reflects a lack of flexibility in the parametrization cho- 
sen in the DSS analysis, albeit not affecting the overall 
quality of the fit. More and better data will be needed 
to constrain more reliably In particular, upcom- 

ing precise data from BELLE [18[ might be sufficient to 
determine the gluon FF from scaling violations in SIA. 
Additional pp data from RHIC and the LHC will add to 
this. We note that as in the case of pions, charm and 
bottom quark FFs and their uncertainties have to taken 
with a grain of salt. 

In order to illustrate the impact of the estimated uncer- 
tainties for the DSS kaon FFs, we compute the charged 
kaon multiplicities at NLO accuracy in the kinematic 
range of the HERMES experiment [24]. To propagate 
the uncertainties of the kaon FFs we use Eq. ([7]) with 
our Hessian sets for both A\ 2 /x 2 = 2% and 5%. The re- 
sults are shown in Fig. [TT] For positively charged kaons, 
the uncertainties resemble in shape those found in the 
case of pions, see Fig. [5j i.e., they are largely indepen- 
dent of z and the energy scale Q. As has to be expected 
from a comparison of Figs. [7J andfTUl the differences be- 
tween results obtained with the Hessian sets for the 2% 
and 5% tolerances are more pronounced for kaons than 
for pions. For negatively charged kaon multiplicities the 
propagated uncertainties behave, however, rather differ- 
ently. We observe an increase in size with the scale Q, 
and, in addition, the estimates for the 2% and 5% Hes- 
sian sets differ much more significantly than for K + and 
7T multiplicities. These features can be associated with 
the larger sensitivity of the K ~ multiplicities to the less 
well determined unfavored kaon FFs as well as with the 
more erratic behavior of the fitted SIDIS data. 

We note that kaon FFs are extremely relevant in con- 
nection with the proper extraction of the strangeness he- 
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FIG. 10: As in Fig. but now showing the uncertainty bands for the NLO DSS kaon FFs at Q 2 = 10 GeV 2 and Q 2 = M|. 




FIG. 11: As in Fig.[S]but for charged kaon multiplicities. On the right-hand-side, only the uncertainty bands for K are visible, 
which are larger than those for K + . 
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licity PDF from spin-dependent SIDIS data [JQ. Our 
error estimates for pion and kaon FFs will be particularly 
useful in future global analyses of the spin structure of 
the nucleon as they allow one to straightforwardly prop- 
agate these sources of uncertainties. Upcoming new data 
on multiplicities for identified hadrons [l9| will be crit- 
ical in further reducing uncertainties of fully flavor and 
charge separated sets of FFs. 

IV. CONCLUSIONS 

We have performed a detailed estimate of the uncer- 
tainties inherent to the extraction of pion and kaon frag- 
mentation functions as performed in the DSS global anal- 
ysis framework. We have carefully examined the validity 
of the approximations underlying the more practical, im- 
proved Hessian method by comparing with the results 
based on the robust Lagrange multiplier technique. 

Even though we have found some differences between 
uncertainties obtained with the two methods, they can 
be readily understood and accounted for by choosing 
are larger tolerance in the Hessian approach. In gen- 
eral, the agreement is much better for the rather well 
constrained parton-to-pion fragmentation functions and 



somewhat less satisfactory for kaons, where data are less 
precise. 

We provide Hessian eigenvector basis sets for both 
pion and kaon fragmentation functions and correspond- 
ing to two different error estimates |23j. These sets 
will greatly facilitate the propagation of uncertainties re- 
lated to fragmentation functions to observables such as 
single-inclusive hadron production cross sections or mul- 
tiplicities. Good knowledge of fragmentation functions 
and their uncertainties is also relevant in understanding 
and analyzing results for spin-dependent semi-inclusive 
deep-inelastic scattering and the modification of identi- 
fied hadron yields in heavy-ion collisions. Our results will 
also prove to be relevant in quantifying the impact of fu- 
ture measurements of pion and kaon yield in upcoming 
global analyses of fragmentation functions. 
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